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Abstract 

We study scaling laws characterizing the inter-diffusive zone between two miscible fluids flowing side by 
side in a Y-shape laminar micromixer using the lattice Boltzmann method. The lattice Boltzmann method 
solves the coupled 3D hydrodynamics and mass transfer equations and incorporates intrinsic features of 3D 
flows related to this problem. We observe the different power law regimes occumng at the center of the 
channel and close to the top/bottom wall. The extent of the inter-diffusive zone scales as square root of the 
axial distance at the center of the channel. At the top/bottom wall, we find an exponent 1/3 at early stages 
of mixing as observed in the experiments of Ismagilov and coworkers [Appl. Phys. Lett. 76, 2376 (2000)]. 
At a larger distance from the entrance, the scaling exponent close to the walls changes to 1/2 [J.-B. Salmon 
et al J. Appl. Phys. 101, 074902 (2007)]. Here, we focus on the effect of finite aspect ratio on diffusive 
broadening. Interestingly, we find the same scaling laws regardless of the channel's aspect ratio. However, 
the point at which the exponent 1/3 characterizing the broadening at the top/bottom wall reverts to the 
normal diffusive behavior downstream strongly depends on the aspect ratio. We propose an interpretation 
of this observation in terms of shear rate at the side walls. A criterion for the range of aspect ratios with 
non-negligible effect on diffusive broadening is also provided. 

PACS numbers: 47.61.Ne, 47.61.Jd, 47.15.Rq, 47.ll.Qr, 87.10.Hk. 
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I. INTRODUCTION 



Microfluidic devices are becoming a means of performing low cost and high-throughput chem- 
ical and biochemical analysis on chip. The range of applications include measuring dynamics 
of protein folding yj], kinetics of enzyme reactions |12|] and surface patterning of cells and pro- 
teins yQ. Some of these applications generally involve cross-stream interaction of two or more 
fluids flowing side by side in a channel. However, such systems are characterized by laminar flow 
due to the small dimensions involved and thus fluids flowing side by side can only mix or interact 
by molecular diffusion. 

In pressure driven flow through rectangular microchannels, where the fluid motion through the 
channel is actuated by pressure pumps, the channel velocity profile is approximately parabolic 
across the shortest dimension, as dictated by the balance between the gradient of the applied pres- 
sure and the viscous shear stress in combination with the no-slip boundary condition at the walls. 
Due to slow fluid motion close to the walls, tracer particles are hardly advected by the flow, thereby 
spending longer time at the walls. At the center of the channel, on the other hand, the flow veloc- 
ity is quite high and tracer particles are transported more efficiently. This results in a fairly wide 
distribution of the amount of time spent by tracer particles in the channel, the so-called residence 
time. As demonstrated experimentally j^]^, this wide distribution of residence times gives rise to 
a non-uniform diffusive broadening across the channel. Ismagilov et al. ||5|], for example, showed 
using confocal fluorescence microscopy that in contrast to the center line of the channel, where the 
extent (along y direction) of the interdiffusion zone exhibits normal diffusive behavior 6 ~ x^^^ 
(x=distance from the inlet), it scales as one third power of the axial distance, S ~ x^^'^, close to the 
top and bottom walls. These observations are shown to be in line with results of a scaling analysis, 
which makes use of a certain similarity with the general Leveque problem [|5|,|6|]. 
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confirms the existence of the 



Numerical solution of the advection-diffusion problem [|7|, 
experimentally observed scaling behavior with an exponent of 1/3 in the proximity of the walls as 
compared to the exponent 1/2 at the center of the channel. Furthermore, these calculations also 
show that the exponent 1/3 can only be observed if the distance (along the flow) from the entrance 
of the channel is not too large. At sufficiently large distances from the inlet, on the other hand, the 
exponent for diffusive broadening close to the walls approaches 1/2 and thus becomes identical 
to the scaling exponent in the center of the channel. The cross over distance, Xc, is identified as 
the distance at which the tracer concentration along the shortest channel dimension (z direction) 
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becomes homogeneous (note that S is measured along the y direction). 

The above mentioned numerical approaches assume a one dimensional parabolic velocity pro- 
file and neglect any dependence of the fluid velocity on the distance x from the inlet as well as on 
the "neutral" direction, y. The first assumption is valid at axial distances x > W, H (W^=width 
and iJ=height of the channel), where the flow is fully developed. The dependence on y, on the 
other hand, can only be neglected if the width of the channel is large compared to the height H, 
i.e. in the case of large aspect ratio, W/H ^ 1. In the experiments |15l|, however, the fully de- 
veloped fluid velocity is two dimensional, i.e. it depends both on y and z due to the finite aspect 
ratios investigated (W/H = 2 — 5) [Qllil. Despite this fact, the agreement between theory and 
experiment is quite good suggesting that the specific form of the velocity profile does not play a 
crucial role as long as a linear regime close to the walls and an approximately uniform flow at the 
channel center can be assumed. 

In the present work, we are going to focus on this aspect via a systematic study of the effect of 
finite aspect ratio on the spreading dynamics of a tracer field entering the Y-junction micromixer 
through one of the arms (see Fig. [I]). For this purpose, we solve, via the lattice Boltzmann (LB) 
method il2|] . the full three dimensional advection-diffusion problem for rectangular channels with 
various aspect ratios. Our previous studies of wall roughness effects on the chaotic mixing of 
passive tracers in a 2D channel showed the flexibility of the LB method in dealing with advection- 
diffusion problem in microchannels IibI [l4 , [isll . The present work extends this approach to 3D 
with a particular focus on various scaling laws within the laminar flow regime using smooth walls 
(no wall roughness effects). 

In addition to a systematic study of the effects related to a finite aspect ratio, our studies differ 
from numerical calculations of Salmon and Ajdari in that we do not make any assumption about 
the shape of the velocity profile. Rather, the velocity profile results from the solution of the Navier- 
Stokes equations for the problem of interest. By doing this, we remain as close as possible to real 
experiments and realize a range of velocity profiles from practically parabolic to those with strong 
deviations from parabolic dependence. 

In the following section, we briefly introduce the simulation scheme and also provide some 
benchmark tests for our LB simulation by comparing our results with known analytical solutions 
for the spreading of a point source in a 3D microchannel both with and without walls. Excellent 
agreement with the analytical solutions is found. We then apply in section |lll] the LB method to 
study the extent of the inter-diffusion zone between two fluids flowing side by side in a Y-shape 
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FIG. 1: Geometry of the micromixer with two inlets at 45°. The origin of the co-ordinate system is at the 
middle of the entrance. Measurements are taken after reaching a fully developed flow 

micromixer. We study the combined effect of Peclet number (defined as the ratio Pe = UH/ D, 
where U is the maximum fluid velocity, H the channel height and D the diffusion coefficient of the 
tracer field) and channel's aspect ratio on the broadening. Our studies include both the upstream 
zone where the non-uniform broadening along the y direction occurs (with an extent S ~ x^^^ close 
to the walls and 6 ~ x^/^ at the center of the channel) and downstream region where diffusion 
has had enough time to homogenize the concentration distribution along the vertical z direction 
whereby leading to homogeneous broadening with a one-half power law both at the walls and in 
the center of the channel. A summary compiles our results. 



II. NUMERICAL METHOD AND ITS VALIDATION 



A. The Lattice Boltzmann Method 

The lattice Boltzmann method [16, 17, [isi [3] is a mesoscopic approach that approximately 
solves the continuum Boltzmann equation in a discrete form by dividing space into a regular 
lattice, time into discrete steps and velocities into a finite number of vectors [|20ll . The number 
and direction of the velocities are chosen such that the resulting lattice is symmetric so as to easily 
reproduce the isotropy of the fluid [|2l|] . The density of the fluid at each lattice site is accounted 
for by a one particle probability distribution fi{r, t), where the subscript i represents one of the 
finite velocity vectors e^, r the lattice site and t the time. During each time step particles stream 
along each velocity vector to a neighboring lattice site and collide locally, conserving mass and 
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momentum in the process. The LB equation describing propagation and collision of the particles 
is given by: 

/,(r + e„ t + 1) - Mr, t) = n.^r, t), (1) 
where fij is the collision operator. 



The most widely used LB method is the lattice BGK model il2|] which approximates the colli- 
sion operator by simplifying it to a single time relaxation towards the local equilibrium distribution 
f^"^. The lattice BGK model is given as: 

/r>,t)-/,(r,t) 



/i(r + e„t + l) -/,(r,t) 



r 



(2) 



where r is relaxation time and the equilibrium distribution f^'^ is closely related to the low Mach 
number expansion of the Maxwellian velocity distribution given as [|22n 

I , . 1 , 1 " 



f-\r, t) = Wip 



1 + — (Ci ■ u) H -(e, ■ 



(3) 



In Eq. dS]), Cs is the sound speed on the lattice and Wi is a set of weights normalized to unity. Near 
equilibrium and in the limit of small Knudsen number (= mean free path/characteristic length 
of problem) the macroscopic Navier Stokes equation can be recovered using Chapman-Enskog 
multiscale analysis MM- The relaxation time r is then found to be related to the kinematic viscosity 

asu = (2r - l)/6. 

The density p, velocity u and pressure p are computed from the distribution function using: 

N N 

P = ^fu P'^ = Y1 f'^''^ P = (4) 

i=0 4=0 

The weights Wi in the equilibrium distribution depend on the number of velocities used for the 
lattice. In this work we use the D3Q15 model which has fifteen velocities with weights Wi given 
as: 

2/9 e, = (0,0,0), i = 0; 
w,,= <l/9 e, = (±1,0,0),(0,±1,0),(0,0,±1) z = L...6; (5) 
1/72 e, = (±l,±l,±l),i = 7....14 
The approach can be extended to simulate diffusing solute by introducing passive tracers into 
the flow field. The tracers are advected with the flow velocity u and follow similar relaxation and 
propagation as the fluid medium. In this study, the tracer field gi does not have any effect on the 
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velocity field (passive tracer limit). The LB equation governing propagation and collision of the 
density (concentration) distribution of the passive tracers is given as: 

gi{r + ei,t + l) - gi{r,t) = , (6) 

Td 

with the equilibrium distribution function, g^'^ir, t) = WiC [1 + (e^ ■ u)/cl]. Here, C is the con- 
centration distribution of the tracers given as C = J2iLo 9i- 

Using Chapman-Enskog multiscale analysis, the advection-diffusion equation can also be re- 
covered in the limit of low Mach number and near equilibrium situation. The relaxation time Td is 
then found to be related to the tracer diffusion coefficient as = {2Td — l)/6. 



B. The Y-shape laminar micromixer 

In the Stokes flow regime, assuming a steady flow with a constant pressure gradient dp/dx 
along the channel, the Navier-Stokes equation can be simplified as dp/dx = riS/'^Ux{x,y, z). 
For a rectangular channel of length L, width W and height H with a fully developed flow 
{ux = Ux{y, z)), this equation can be solved via the method of separation of variables with no- 
slip boundary condition on the walls of the channel. One thus obtains the following result for the 



velocity field UM 



dp/dx 
2r]L 



i/2 



— z' 



8H^ 



- 1 cosh(^) /nn 

22 — — - sill 



cosh {^) V H 



z + 



H 



(7) 



n=l,3,5 

Far away from the side walls and close to the center of the channel (y ^ 0), the fluid velocity 
follows the well known 2D parabolic profile. In the Y-shape micromixer studied here, the velocity 
profile at the two inlets is quite different from that in the mixing channel. Assuming that the 
channel width W is larger than its height, W > H,it takes approximately a distance of the order 
of W from the junction for velocity profile to become fully developed. 

In the steady state, the mass transport of a solute of concentration C with a constant isotropic 
diffusion coefficient D in a unidirectional and one dimensional velocity field, u = {ux{z), 0, 0), 
can be describe by 



Ux{z) 



dC_ 

dx 



D 



d^ d^ ^ 

Qrj,2 Qy2 Q^2 



(8) 
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Given that the initial concentration at the junction, where the two fluids meet, is C(x = 0) = Cq 
and the characteristic fluid velocity is U, Eq. ([8]) can be put into a non-dimensional form using 
relevant characteristic length scale. A convenient form is to use x* = x/ (PeH), y* = y/H and 
z* = z/H. Equation ([8]) then becomes 

The first term on the right hand side of Eq. ^ represents the contribution of axial diffu- 
sion. It can be neglected if two conditions are fulfilled. The first condition is that the first term 
on the right hand side must be small compared to the other two terms i.e. d'^C/dx*'^ /Fe^ <^ 
d'^C/dy*^, d'^C/dz*'^. Obviously, this happens when Pe ^ 1. The second condition is that the 
axial convective term must be large compared to the axial diffusive term, i.e. Ux{z*)dC / dx* 3> 
{d'^C/dx*'^) /Pe^. In this case, using order of rnagnitude analysis, the regimes (x*, z*) where this 
condition is valid can easily be identified to be [SI], <^ 1/ {Fe^Ux{z*)). 

Assuming the validity of both the conditions stated above, Eq. ^ can be simplified to 



Salmon and Ajdari [|8|] performed a numerical integration of Eq. (flOl) assuming a parabolic 
velocity profile across the channel, i.e. = Ux{0)[l — {2z/HY]. The lattice Boltzmann method, 
on the other hand, solves the full Eq. ([8]) incorporating all relevant 3D features such as the fact 
that, in general, the fluid velocity is not parabolic and depends on all the three coordinates x, y and 



z, i.e. Ux = Uxix,y,z) MM- 



C. Simulation details 

All simulations reported in this work are conducted using the D3Q15 LB model on a parallel 
machine. The lattice Boltzmann code is written in c-i-i- with MPI (message passing interface) 
routines that enable exchange of information between processes. The code is run on SMP 8 x 
2 core opteron 2.4 GHz. Typical channel dimensions used vary in the range 40 x 20 x 1000 to 
400 X 20 X 1000 as well as a channel of size 40 x 40 x 1000 (lattice units), thereby allowing the 



study of aspect ratios between 1 and 20. At the walls of the channel, the bounce-back scheme [|19|] 
is implemented. Populations streamed to the walls are simply reversed back along the directions 
where they came from. Starting with the fluid at rest, a flow is imposed by the application of a 



body force. Note that it takes a finite amount of time (of the order of the momentum diffusion time 
tdifr = H"^ / (8z/) until steady state is reached with respect to the fluid velocity, i.e. until the 
velocity field becomes independent of time. In order to avoid this transient effects, we wait a time 
of 5 X tdiff before injecting the tracer field. 

After a stationary flow has been reached, we start to continuously inject passive tracers into the 
channel. This is simply achieved by setting at each time step the concentration of the tracers at the 
inlet to 1 (i.e. C(x = 0; t) = Co = 1 for all times t). We then monitor the time evolution of the 
tracer concentration field along the channel. Computation of the extent of diffusive broadening 
of the concentration distribution at a given cross section is done after a time independent 
concentration profile has been reached at the area of interest (corresponding to a dynamic balance 
between the incoming and outgoing tracer populations). The modeled diffusing analyte in this 
work has a diffusion coefficient of £) = 0.5 x 10~^m^/s corresponding to At, where we 

have chosen the lattice units of Ax = 3yum and At = 600ns. The kinematic viscosity of the fluid 
was set to z/ = O.lAx^/At corresponding to the viscosity of water (z/watcr ~ lO^^m^/s at room 
temperature). The parameters chosen are readily comparable with that of the experiments. 



D. A point source in a 3D channel 

Before presenting in section UlI] results of our simulations on laminar micromixer, we provide 
here a simple test of our lattice Boltzmann approach. For this purpose, we compute the funda- 
mental problem of diffusion of a point source placed at point vq = (xq, yo, Zq) at time t = on a 
3D square lattice using the LBGK model. Assuming a constant and isotropic diffusion coefficient 
D, the system is described by the well known diffusion equation dC/dt = DAC (A=the Laplace 
operator) with the initial condition C(r,0) = 6{x — Xo)S{y — yn)5{z — zq). The fundamental 



solution or Green function of this equation is well known to be [|25l] 



In the presence of non-absorbing impermeable walls (equivalent to no-flux boundary condi- 
tion), the analytical solution for the diffusion of a point source can be obtained by using the princi- 
ple of superposition and method of images [2^. The positions (x^, ?/m, Zn) of the infinite number 
of images formed on each side of the box is obtained from the real point source according to the 
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FIG. 2: Comparison of lattice Boltzmann simulation with analytical results (a) diffusion of a point source 
within a cubic box surrounded by non-absorbing walls (equivalent to a no flux boundary condition), (b) 
diffusion of a single point source in a fluid moving with a uniform velocity. The dashed line in (a) represents 
the solution of the same problem in an infinite system (no walls) thus underlying the non-trivial effect of 
the walls, fairly well captured by the LB method. 

relations Xk = Xq + kL^, Vm = yo + rriLy and Zn = Zq + nL^. The final result for the concentration 
profile in the presence of non-absorbing walls is thus 

k,m,n=—QO 

Using a square lattice of dimension Lx=SO, Ly=SO, Lz=SO and diffusion coefficient D = 0.5 



(lattice units), we perform lattice Boltzmann simulation with a simple bounce back scheme Ill9f] 
at the walls. Figure [2ta) compares the simulation results for the concentration profile along the 
x-direction with the analytical solution obtained from Eq. (fT2l) . In the same figure we also plot 
the Gaussian function given by Eq. (fTTI) . which represents the solution of the same problem in the 
absence of the walls, thus emphasizing the non-triviality of the obtained solution in the presence 
of the walls. Our simulation results are in agreement with the analytical results obtained from Eq. 
(fT2l) within an error of 0.2%. This agreement shows that LB model is able to capture the effect of 
the wall correctly. It is straight forward to extend the situation to the case of diffusion in a uniform 
steady unidirectional velocity field u = (Uq, 0, 0) . The results for this situation are shown in 
Fig. |2b) for the case of an infinite system (no walls). Here, the point source exhibits a normal 
diffusive broadening {S ~ t^^^), while at the same time being advected by the flow. The analytic 
curves shown in Fig. |2l;b) are flow advected version of Gaussian distributions given in Eq. (fTTI) . 



where x — xq is replaced by x — (xq + Uot). Noting that the center of the Gaussian distribution 
along the flow direction obeys (x) = Uot, the normal diffusive broadening can also be expressed 



ni. RESULTS AND DISCUSSION 

For a Peclet number of Pe = 1000, we plot in Fig.[3l the cross-sectional image of the concentra- 
tion field at various axial positions x* = x/ (if Pe) along the channel for aspect ratios of W/ H = 2 
and 5 respectively. First note that the width of the interdiffusion zone increases with increasing 
axial distance (compare panels (a), (b) and (c)). This is reminiscent of the role played by time in 
diffusion process. Indeed, recalling that the tracer concentration at a height z is advected with the 
velocity u{z) and neglecting the effect of shear for the moment, a rough correspondence between 
diffusion time and the distance x from the inlet can be obtained via t(x, z) = x/u{z). A larger 
distance from the inlet thus corresponds to a larger diffusion time. 

A survey of Fig. [3] allows a second important observation, namely that the extent 5 of interdif- 
fusion zone increases when going from the center of the channel towards the walls, the so called 
"butterfly effect" (see e.g. panel (b)). Again, making use of the above estimate of diffusion time 
as t(x, z) = x/u{z), a qualitative understanding of this behavior can be gained by noting that u{z) 
is maximum at the center of the channel (z = 0) and decreases to zero at the walls {z = ±H/2). 
Thus, at a given distance x from the inlet, the time available for diffusion is larger in the proximity 
of the walls as compared to the center of the channel. 

Even though being able to describe some important qualitative features of diffusive broadening, 
the above argument is too crude to capture the different scaling laws discussed above (6 oc x^/^ 
close to the walls as compared to 5 oc x^^^ at the channel center). Indeed, an argument based only 
on an estimate of the effective diffusion time would yield 6 = VODt = y^6Dx/u{z) oc x^/^. 
In other words, one would obtain diffusive broadening with an exponent of 1/2 for all distances 
from the wall, the z-dependence being reflected in an effective diffusion coefficient D/u{z). An 
adequate description of experimental observation, therefore, requires taking into account the effect 
of the non-uniformity of the flow on tracer distribution. Interestingly, assuming a linear velocity 
profile seems to be sufficient for a derivation of the observed scaling exponent of 1/3 close to the 
walls isi Isl]. The exponent 1/2, on the other hand, results from the presence of a quasi-uniform 
flow in the central region of the channel. 
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A third important aspect in Fig. [3] is the effect of aspect ratio. A comparison of the left and 
right panels of Fig.[3]at a given distance from the inlet clearly shows that the interdiffusion zone is 
broader in the case of lower aspect ratio. As a consequence, also the inhomogeneity of the tracer 
concentration field along the z direction is more pronounced when the aspect ratio is decreased. 
One could therefore expect that a longer diffusion time is required in order to homogenize the 
concentration field across the vertical direction. Using the rough correspondence between time and 
axial distance along the lines discussed above, the position at which the concentration distribution 
across z direction becomes homogeneous should increase for smaller aspect ratios. A comparison 
of panels (c) and (f) in Fig. [3] confirms this expectation. Both these panels correspond to the same 
distance from the inlet but different aspect ratios. While in Fig. [Sj^f) (aspect ratio=5) the solute 
concentration is already homogeneous along the ^-direction, it exhibits significant inhomogeneity 
in Fig.[3];c) (aspect ratio=2). An important consequence of this feature will be worked out below. 

In order to proceed with a more quantitative analysis, we adopt a simple definition of the 
extent of interdiffusion zone b as the width of the concentration profile at which the concen- 
tration is reduced to 20% of the maximum value at the inlet. It is to be stressed that the re- 
sults presented here are insensitive to the specific definition of b. Other definitions such as 
using the integral over the second moment of the spatial derivative of the concentration field, 
S = (J dyy'^dC / dy) / (J dydC/dy), lead essentially to the same conclusions. The present sim- 
ple definition, however, is numerically more robust since it does not require the computation of 
numerical derivative of the concentration field at discrete intervals. 

As discussed in section UIl 5 ~ x^^"^ for diffusion in a uniform velocity field. In the case 
of the inhomogeneous velocity profile considered here, the relation between the extent of the 
broadening 6 and the axial distance x takes the more general form 6 x'^, where 7 is the exponent 
characterizing the diffusive broadening. It is shown in Ref. Js] that 7 in general depends both on 
X and z. In order to examine this property within our simulations, we survey in Fig.|4]the change 
of 6 downstream along the x direction. The local slope of the (log-log) plots gives the scaling 
exponent 7. 

Within each panel in Fig.lH the extent S of the interdiffusion zone is shown versus axial distance 
both for the center of the channel (z/H = 0) as well as in the proximity of one of the walls 
(z/H = 0.45). The following features can be observed in the two left panels shown in Fig. |4] 
corresponding to a Peclet number of Pe = 1000. At small distances from the inlet, 5 ~ x^/^ at the 
center of the channel while 6 ~ x^^'^ close to the wall. However, increasing the lattice resolution 
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FIG. 3: Solute concentration C iny — z cross sections of a rectangular channel for aspect ratios W/H = 2 
(left panels) and 5 (right panels) at various reduced distances x* from the inlet. From top to bottom: 
X* = 10~^'®, 10~^'^ and 10~^'^^. For all cases shown, the Peclet number is Pe = 1000 and the height of 
the channel is i7 = 20 lattice units. The width of the channel is 1^ = 40 (left panels) and W = 100 (right 
panels) lattice units. 
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FIG. 4: Plot of log{6/H) versus log(a;*) for aspect ratios of 2 (upper panels) and 5 (lower panels) at Peclet 
numbers of Pe = 1000 (left panels) and Pe = 10000 (right panels). The local slope of the lines gives the 
scaling exponents. The .x-axis is non-dimensionalized via x* = x/{HPe). The channel height is = 20 
(left panels) and H = 200 (right panels). The channel widths ai^e chosen such that the aspect ratio is 
W/ H = 2 in the case of upper panels and W/ H = 5 in the case of lower panels (W = 40 (top left), 
W = 400 (top right), W = 100 (bottom left) and W = 1000 (bottom right)). In all the panels shown, the 
vertical dashed line indicates x* = 1/8 as obtained from a dimensional estimate of the distance for vertical 
homogenization of the solute concentration. 

of our simulation domain in the z direction and correspondingly the Peclet number reveals another 
short regime, where 5 ~ x^/^ across the entire cross section of the channel. This regime, as shown 
in Fig.|4]for Peclet number of Pe = 10000, occurs earlier in the flow at the entrance of the channel 
and changes to the 1/3 regime over a short time interval. 

The transition of this 1/2 regime to the 1/3 regime can be understood using an analytical argu- 
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(a) (b) (c) 




FIG. 5: Plot of the scaling exponent 7 (appearing in 5 oc x*'^) in the log(x*) — z* plane for aspect ratios of 



(a) 5 and (b) 2. In (c) the difference, 7(b)-7(a), is shown. The exponent 7 is obtained by local fits to log-log 
plot in Fig. m 




Aspect ratio (W/H) Aspect ratio (W/H) 

FIG. 6: Left: Plot of the cross over point Xc versus the aspect ratio for different Peclet numbers as indicated. 
Right: The same data as in the left panel, rescaled via x* = Xc/{HPe). The horizontal solid hne marks 
the value of x* = 1/8 obtained from an estimate of the time necessary for vertical homogenization in a 
channel with infinite aspect ratio (W/H —>■ 00): 2Dtc = Using the cross-over distance Xc and the 

mid-channel fluid velocity U to estimate the cross-over time, = Xc/U, one obtains Xc = H^U/{8D) 
and hence x* = 1/8 (recall our definition of the Peclet number Pe = HU/D). This estimate is improved 
significantly by taking, tc = Xc/U, where U is the average fluid velocity across the channel (dashed 
horizontal hne). 
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FIG. 7: (a) Velocity, Ux, and shear rate, Fz = dux/dz, versus y for aspect ratios of W/H = 2 and 5. The 
distance from the bottom wall is equal to z = H/10 (=2 lattice units) for all the data shown. For the aspect 
ratio of 5, there is a wide range of y-values for which both Ux and are roughly constant thus justifying 
the assumption Ux{y, z) ~ Ux{z). This assumption, however, fails at the smaller aspect ratio shown. In 
this case, there is practically no y independent region. Vertical dashed lines mark y = and y = SH/A 
for which velocity profiles along the vertical (z) direction are depicted in the adjacent panel, (b) Velocity 
versus z for y = and y = 3H/4 for aspect ratios of W/H = 2 and 5. In the case of an aspect ratio of 
5, the velocity profile is identical for both values of y. This is in accordance with the panel (a) where Ux 
hardly varies in the y range delimited by the two vertical dashed lines. At a lower aspect ratio of 2, on the 
other hand, the effect of the side wall is quite significant. The Peclet number is Pe = 1000 in all the cases 
shown. 

ment similar to Leveque analysis by assuming a linear velocity profile very close to the top/bottom 
wall Isl. Qualitatively, this cross-over represents the enhancement in diffusive broadening arising 
from a homogeneous shear rate. Since, in the case of the Poiseuille flow studied here, the shear 
rate is practically zero in the center of the channel, the effect appears only close to the walls, where 
the velocity profile is approximately linear. 

At larger distances from the inlet, on the other hand, the scaling exponent close to the wall 
gradually increases towards the value of the exponent at the center of the channel, the latter being 
very close to 1/2. This latter behavior can be understood by assuming that the cross over from 
an exponent of 7 = 1/3 to 7 = 1/2 corresponds to a homogeneous tracer distribution along the 
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2;-axis. A criterion for vertical homogenization via diffusion is obtained from 2Dtc = {H/2Y 
which, using the maximum fluid velocity U to estimate the cross over time, tc = Xc/U, yields the 
cross over distance Xc = H'^U/{8D) = HFe/8. In terms of the reduced distance, this relation 
translates to x* = 1/8 [SD. 

The above estimate of Xc does not take into account the effect of aspect ratio. In fact, as 
discussed above, the vertical homogenization takes place at larger axial distances when the aspect 
ratio is decreased (compare panels (f) and (c) in Fig. [3]). In order to investigate this issue, data in 
Fig.|4]are plotted for two different aspect ratios of W/H = 2 (upper panels) and W/H = 5 (lower 
panels). Indeed, a comparison of the panels (a) and (b) [as well as (c) and (d)] in Fig. |4] suggests 
that Xc increases when the aspect ratio is decreased. The left and right panels in Fig. |4] differ in 
Peclet number investigated. This is to underline the fact that the observed trend with regard to 
aspect ratio is not related to the specific choice of the Peclet number. 

In order to quantify the effect of aspect ratio further, we determine for two aspect ratios of 
W/ H = 2 and W/ H = 5 the values of the exponent 7 for each position x and z along the channel 
by performing a running fit on the curves in Fig.|4]as proposed in [18|]. The results are shown in 
Fig. [51 The dotted black lines in Fig.|5l^) indicate concentration boundary layer which grows with 
the axial distance at a rate z* ~ x*^^^ The growth of the concentration boundary layer results 
from the diffusive flux of solute from the top/bottom wall to the center of the channel. The white 
area of the plot corresponds to the high exponent region of the channel. Our resolution is limited 
by the small number of points in the z* direction. The difference between the exponents obtained 
for the two investigated aspect ratios is shown in Fig. [3c) confirming the retarded cross over to 
normal diffusive behavior in the case of smaller aspect ratio. 

In a more systematic study, we varied the aspect ratio from 1 to 20 and determined the cross over 
distance Xc- The left panel of Fig. [6] illustrates the behavior of the thus obtained cross over point 
within the range of aspect ratios investigated and for different Peclet numbers of Pe = 500, 1000 
and 10000. In line with the results presented above, we observe an increase in the cross over point 
as the aspect ratio decreases. The right panel of Fig. [6] depicts a non-dimensionalized version of 
the same data. Interestingly, the data for different Peclet numbers collapse onto a master curve in 
the limit of high aspect ratios. This is an important observation, since it suggests that, in this limit, 
the effect of Peclet number is indeed a mere rescaling of time or axial distance. 

In an attempt to better understand the reason for dependence of Xc on aspect ratio, we plot 
in Fig. Ulsi) fluid velocity as well as its spatial derivative, = du^/dz, as a function of the 
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transverse coordinate y for aspect ratios of W/H = 2 and W/H = 5. The plot is done for a 
distance of 2 lattice units from the bottom wall. For the aspect ratio of 5, there is a wide range 
of y-values for which both and are roughly constant thus justifying the assumption that 
u,j.{y, z) = Ux{z) (independent of y). This assumption, however, fails at the smaller aspect ratio 
shown. The data shown in the right panel of Fig. |7] further underline the importance of taking into 
account the dependence of the fluid velocity both on z and on y when the aspect ratio is small. 

The data shown in Fig. |7] suggest that the increase in Xc is probably due to the side wall shear 
effect on solute distribution. At a lower aspect ratio, the shear rate decreases strongly close to 
the side wall. An estimate from Fig. Ul^) shows that for an aspect ratio of 2, the at position 
y = 3H/4: is about 40% less than that at the center of the channel (y = 0). Therefore, solute 
diffusing towards the side walls samples a strongly decreasing shear rate. Given that the extent of 
broadening 6 at the top/bottom wall takes the form 5 ~ {xD /F^) |5l] , the interdiffusion becomes 
faster when approaching the side walls, whereby enhancing the inhomogeneity of diffusion (the 
so called "butterfly effect") which is already present at infinite aspect ratio. This is exactly what 
we observe in the concentration profile images shown in Fig.O 

In general, solute diffusing towards the side walls from the center of the channel samples con- 
stant velocity gradient up to a distance H from the side wall. The time spent by solute before 
sampling a substantial decrease in the velocity gradient is, therefore, tw ~ {W/2 — /2D. For 
W ^ H, this time scale is greater than the time scale to diffuse from the top/bottom wall to the 
center of the channel denoted as tn ~ {H/2Y /2D. Thus, at high aspect ratio, solute concentration 
becomes homogeneous along the vertical direction long before the side walls are "felt". Conse- 
quently, the cross over point Xc becomes independent of aspect ratio. In the case where tw < ^h, 
on the other hand, one can not neglect the additional enhancement of inhomogeneity of diffusive 
broadening due to the side wall effect. This leads to < 3i/ for the effect of aspect ratio being 
significant. As a survey of x* in Fig. [6tb) reveals, this simple estimate (which is based on a di- 
mensional argument only) lies within a factor of two of the result obtained within our computer 
simulations. 

A. Summary 

We study the effect of a finite aspect ratio on transverse diffusive transport of miscible so- 
lutes flowing in a pressure driven microchannel using the lattice Boltzmann method. The lattice 
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Boltzmann method incorporates the essential 3D features such as the non-parabolicity of the ve- 
locity profile and the velocity gradient due to the side walls. We observe the previously reported 



n 



n 



ISi LZL ISL |9[| different exponents characterizing the extent of the broadening both at the early stage 
of mixing of the two fluids and at the later stage downstream. Interestingly, we observe the same 
scaling laws regardless of the channel aspect ratio. However, extent of diffusive broadening and 
the position, Xc, at which the broadening becomes uniform and finally reverts to the 1/2 behavior 
vary remarkably with the channel aspect ratio. The Peclet number, on the other hand, is found to 
play the role of a scale factor in Xc in a way that x* = Xc/ (HFe) is independent of aspect ratio. 
This is inline with the general structure of the advection-diffusion equation, upon neglection of 
axial diffusion. 

A qualitative understanding of the effect of aspect ratio is provided invoking the influence of 
the shear stress in the proximity of the side walls on diffusive broadening. This is based on the 
idea that side wall shear stress is non-negligible in a region of width H close to the side walls. 
The corresponding effect on inhomogeneous diffusive broadening will be felt if the solutes at the 
center of the channel have enough time to reach this region before vertical homogenization takes 
place. This allows to derive a simple criterion to decide whether a given aspect ratio is "large" or 
"small". Within a factor of two, this estimate correctly reproduces the behavior observed within 
our lattice Boltzmann computer simulations. 
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